# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import itertools as it
import numpy as np
import scipy as sp
import sympy as sm
import warnings
try:
import flint
has_flint = True
except ImportError:
import warnings
from hysop.tools.warning import HysopPerformanceWarning
msg = "Failed to import python-flint module, falling back to slow sympy solver."
warnings.warn(msg, HysopPerformanceWarning)
flint = None
has_flint = False
from hysop.tools.henum import EnumFactory
from hysop.tools.htypes import check_instance, InstanceOf, to_tuple
from hysop.tools.cache import update_cache, load_data_from_cache
from hysop.tools.io_utils import IO
from hysop.tools.sympy_utils import tensor_xreplace, tensor_symbol
from hysop.tools.decorators import debug
from hysop.tools.warning import HysopCacheWarning
from hysop.numerics.stencil.stencil_generator import CenteredStencilGenerator, MPQ
from hysop.numerics.interpolation.interpolation import MultiScaleInterpolation
def _check_matrices(*x):
return all(_check_matrix(xi) for xi in x)
def _check_matrix(x):
return np.isfinite(np.asarray(x).astype(np.float64)).all()
PolynomialInterpolation = EnumFactory.create(
"PolynomialInterpolation",
[
"LINEAR", # requires 0 ghosts (no derivatives required)
"CUBIC", # derivatives order is specified with SpaceDiscretization
"QUINTIC", # derivatives order is specified with SpaceDiscretization
"SEPTIC", # derivatives order is specified with SpaceDiscretization
"NONIC", # derivatives order is specified with SpaceDiscretization
"CUBIC_FDC2", # requires 1 ghosts (estimate derivatives with 2nd order centered fd)
"CUBIC_FDC4", # requires 2 ghosts (estimate derivatives with 4th order centered fd)
"CUBIC_FDC6", # requires 3 ghosts (estimate derivatives with 6th order centered fd)
"QUINTIC_FDC2", # requires 1 ghosts
"QUINTIC_FDC4", # requires 2 ghosts
"QUINTIC_FDC6", # requires 3 ghosts
"SEPTIC_FDC2", # requires 2 ghosts
"SEPTIC_FDC4", # requires 3 ghosts
"SEPTIC_FDC6", # requires 4 ghosts
"NONIC_FDC2", # requires 2 ghosts
"NONIC_FDC4", # requires 3 ghosts
"NONIC_FDC6", # requires 4 ghosts
],
)
[docs]
class PolynomialInterpolator:
[docs]
@classmethod
def build_interpolator(cls, pi, dim, fd=None, verbose=False, approximative=False):
check_instance(pi, PolynomialInterpolation)
kwds = {"dim": dim, "verbose": verbose, "approximative": approximative}
spi = str(pi)
if spi.startswith("LINEAR"):
deg = 1
fd = 2
elif spi.startswith("CUBIC"):
deg = 3
elif spi.startswith("QUINTIC"):
deg = 5
elif spi.startswith("SEPTIC"):
deg = 7
elif spi.startswith("NONIC"):
deg = 9
else:
msg = f"Unknown PolynomialInterpolation value {pi}."
raise NotImplementedError(msg)
for i in range(1, 5):
if spi.endswith(str(2 * i)):
fd = 2 * i
break
if fd is None:
msg = "Could not determine finite differences order."
raise RuntimeError(msg)
obj = cls(deg=deg, fd=fd, **kwds)
obj.spi = spi
return obj
[docs]
@classmethod
def cache_file(cls):
_cache_dir = IO.cache_path() + "/numerics"
_cache_file = _cache_dir + "/polynomial.pklz"
return _cache_file
def __init__(self, dim, deg, fd, approximative=False, verbose=False):
"""
Create a PolynomialInterpolator.
Parameters
----------
dim: int
Number of dimensions to interpolate.
deg: int or tuple of ints
Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
Degree should be odd on each axis.
fd: int or tuple of ints
Order of centered finite differences used to compute derivatives in each direction.
Will affect the number of ghosts of the method.
Should be even because this interpolator only use centered dinite differences.
approximative: bool
Use np.float64 instead of exact fractions to compute weights.
verbose: bool
Enable or disabled verbosity, default to False.
Attributes
----------
dim: int
Number of dimensions to interpolate.
deg: tuple of ints
Polynomial degree (1=linear, 3=cubic, 5=quintic, 7=septic, ...)
Degree should be odd: deg=2k+1
fd: tuple of ints
Order of centered finite differences stencils used to compute derivatives for
each direction.
p: tuple of ints
Corresponds to deg+1.
The total number of polynomial coefficients corresponds to P=p0*p1*...*p(dim-1).
P: int
The total number of polynomial coefficients P=p0*p1*...*p(dim-1)
k: tuple of ints
Max derivative order required to compute the polynomial interpolator coefficients
in each direction.
Also the regularity of the resulting interpolant. Corresponds to (deg-1)/2.
ghosts: tuple of ints
Return the number of ghosts required by the interpolator on each axis.
Corresponds to (k>0)*[fd//2 - 1 + (k+1)//2]
deg k (k+1)/2 | FDC2 FDC4 FDC6
linear: 1 0 0 | 0 0 0
cubic: 3 1 1 | 1 2 3
quintic: 5 2 1 | 1 2 3
septic: 7 3 2 | 2 3 4
nonic: 9 4 2 | 2 3 4
n: tuple of ints
Corresponds to 2*(ghosts+1), the number of required nodes to generate the
polynomial coefficients (in each direction).
In total we have N=n0*n1*...*n(dim-1) input nodes.
G1 G1
<-> <->
X X X X
X P P X
X P P X
X X X X
<----->
n1
N: int
Total number of input nodes N=n0*n1*...*n(dim-1).
M: np.ndarray
Grid values to polynomial coefficient matrix:
M.dot(F.ravel()) will give C.ravel(), coefficients of P(x0,x1,...)
N
<--------->
X X X X X X ^ |f0| ^ |c0| ^
X X X X X X | |f1| | |c1| |
M = X X X X X X | P F = |f2| | N C = M*F = |c2| | P
X X X X X X | |f3| | |c3| |
X X X X X X v |f4| | |c4| v
|f5| v
If approximative is set to True, M will contain np.float64
Else is will contain rationals.
See Also
--------
:class:`PolynomialSubgridInterpolator`: Precompute weights for fixed subgrid
interpolation.
"""
assert dim > 0, "dim<=0"
deg = to_tuple(deg)
if len(deg) == 1:
deg *= dim
check_instance(deg, tuple, values=int, size=dim)
fd = to_tuple(fd)
if len(fd) == 1:
fd *= dim
check_instance(fd, tuple, values=int, size=dim)
p = tuple(degi + 1 for degi in deg)
k = tuple((degi - 1) // 2 for degi in deg)
ghosts = ()
n = ()
for fdi, ki in zip(fd, k):
if ki > 0:
gi = (fdi // 2) - 1 + (ki + 1) // 2
else:
gi = 0
ni = 2 * (gi + 1)
ghosts += (gi,)
n += (ni,)
check_instance(deg, tuple, values=int, size=dim)
check_instance(fd, tuple, values=int, size=dim)
check_instance(p, tuple, values=int, size=dim)
check_instance(k, tuple, values=int, size=dim)
check_instance(n, tuple, values=int, size=dim)
check_instance(ghosts, tuple, values=int, size=dim)
assert all(fdi % 2 == 0 for fdi in fd), "fd % 2 != 0"
assert all(degi % 2 == 1 for degi in deg), "deg % 2 != 1"
assert all(pi % 2 == 0 for pi in p), "p % 2 != 0"
assert all(ni % 2 == 0 for ni in n), "n % 2 != 0"
P = np.prod(p, dtype=np.int32)
N = np.prod(n, dtype=np.int32)
self.dim = dim
self.deg = deg
self.p = p
self.P = P
self.k = k
self.n = n
self.N = N
self.fd = fd
self.ghosts = ghosts
self.approximative = approximative
self.verbose = verbose
self.key = ("PolynomialInterpolator", dim, deg, fd, approximative)
self._build_interpolator()
def _collect_stencils(self):
dim = self.dim
ghosts = self.ghosts
verbose = self.verbose
fd = self.fd
approximative = self.approximative
k = self.k
n = self.n
ghosts = self.ghosts
if verbose:
print("\nCollecting 1D stencils:")
SG = CenteredStencilGenerator()
SG.configure(dim=1, dtype=np.float64)
S = {}
for direction in range(dim):
if verbose:
print(f" Direction {direction}")
Sd = S.setdefault(direction, [])
nd = n[direction]
kd = k[direction]
fdd = fd[direction]
gd = ghosts[direction]
for i in range(kd + 1):
msg = "Failed to compute stencil derivative={}, order={}, origin={}"
msg = msg.format(i, fdd, gd)
try:
if approximative:
Si = SG.generate_approximative_stencil(order=fdd, derivative=i)
else:
Si = SG.generate_exact_stencil(order=fdd, derivative=i)
Si.replace_symbols({Si.dx: 1})
except:
print(msg)
raise
msg += f" got {Si.coeffs}."
assert not Si.is_symbolic(), msg
Si = Si.reshape((nd - 1,))
assert Si.origin == gd
Si = Si.coeffs
Sd.append(Si)
if verbose:
print(f" {i}-th derivative: {Si}")
return S
def _build_stencil(self, dvec):
dvec = np.asarray(dvec)
k = self.k
S = self.S
assert dvec.size == self.dim, "dvec.size != dim"
assert all(dvec >= 0), f"dvec < 0 => {dvec}"
assert all(di <= ki for (di, ki) in zip(dvec, k)), f"dvec > dmax => {dvec}"
Sd = S[0][dvec[0]].copy()
for d, i in enumerate(dvec[1:], 1):
Sd = np.tensordot(Sd, S[d][i], axes=0)
return Sd
def _build_interpolator(self):
dim = self.dim
deg = self.deg
k = self.k
n = self.n
p = self.p
ghosts = self.ghosts
approximative = self.approximative
verbose = self.verbose
xvals, xvars = tensor_symbol("x", shape=(dim,))
fvals, fvars = tensor_symbol("F", n, ghosts)
pvals, pvars = tensor_symbol("C", p)
self.xvals, self.xvars = xvals, xvars
self.fvals, self.fvars = fvals, fvars
self.pvals, self.pvars = pvals, pvars
try:
data = load_data_from_cache(self.cache_file(), self.key)
if data is not None:
(P0, S, M) = data
if _check_matrix(M):
self.P0 = P0
self.S = S
self.M = M
return
except Exception as e:
msg = f"Failed to load data from cache because:\n{e}"
warnings.warn(msg, HysopCacheWarning)
P0 = 0
for idx in it.product(*tuple(range(0, pi) for pi in p)):
P0 += pvals[idx] * np.prod(np.power(xvals, idx))
self.P0 = P0
S = self._collect_stencils()
self.S = S
if verbose:
print("\nGenerating variables:")
print(" *space vars: ")
print(xvals)
print(" *grid values:")
print(fvals)
print(" *polynomial coefficients:")
print(pvals)
print(" *polynomial patch:")
print(P0)
print("\nBuilding system...")
eqs = []
for dvec in it.product(*tuple(range(0, ki + 1) for ki in k)):
if verbose:
print(f" => derivative {dvec}")
dP0 = P0
for i, deg in enumerate(dvec):
dP0 = sm.diff(dP0, xvals[i], deg)
stencil = self._build_stencil(dvec)
if verbose:
print(" stencil:")
print(stencil)
for idx in it.product(*tuple(range(gi, gi + 2) for gi in ghosts)):
if verbose:
print(f" -> point {idx}")
pos = np.asarray(idx) - ghosts
pos = dict(zip(xvals, pos))
eq = dP0.xreplace(pos)
for offset in it.product(*tuple(range(-gi, gi + 1) for gi in ghosts)):
fidx = tuple(np.add(idx, offset))
sidx = tuple(np.add(offset, ghosts))
eq -= fvals[fidx] * stencil[sidx]
eqs.append(eq)
if verbose:
print(f" {eq}")
# Build system such that A*c = B*f where c are the polynomial coefficients and
# f the node values
dtype = np.float64 if approximative else object
A = np.empty((self.P, self.P), dtype=dtype)
B = np.empty((self.P, self.N), dtype=dtype)
assert len(eqs) == self.P
for i, eq in enumerate(eqs):
for j, ci in enumerate(pvars):
A[i, j] = +eq.coeff(ci)
for j, fi in enumerate(fvars):
B[i, j] = -eq.coeff(fi)
# C = Ainv*B*f = M*f
if verbose:
print("\nSolving system...")
if approximative:
Ainv = np.linalg.inv(A)
elif has_flint:
coeffs = list(flint.fmpq(x.p, x.q) for x in A.ravel())
Afmpq = flint.fmpq_mat(*(A.shape + (coeffs,)))
Afmpq_inv = Afmpq.inv()
coeffs = list(
sm.Rational(sm.Integer(x.p), sm.Integer(x.q))
for x in Afmpq_inv.entries()
)
Ainv = np.asarray(coeffs).reshape(A.shape)
else:
# /!\ sympy is really slow
Ainv = np.asarray(sm.Matrix(A).inv())
if verbose:
print("\nBuilding matrix...")
M = Ainv.dot(B)
self.M = M
update_cache(self.cache_file(), self.key, (P0, S, M))
[docs]
def interpolate(self, fvals):
"""Return the polynomial interpolating input node values"""
fvals = np.asarray(fvals)
assert fvals.shape == self.fvals.shape
pvals = self.M.dot(fvals.ravel())
P0 = self.P0.xreplace(dict(zip(self.pvars, pvals)))
return sm.utilities.lambdify(self.xvars, P0)
[docs]
def generate_subgrid_interpolator(self, grid_ratio, dtype=None):
return PolynomialSubgridInterpolator(
interpolator=self, grid_ratio=grid_ratio, dtype=dtype
)
def __hash__(self):
objs = (self.dim, self.deg, self.fd, self.approximative)
return hash(objs)
[docs]
class PolynomialSubgridInterpolator:
def __init__(self, interpolator, grid_ratio, dtype=None):
"""
Create a PolynomialSubgridInterpolator from a PolynomialInterpolator and a number of
subrid points.
Parameters
----------
interpolator: PolynomialInterpolator
Interpolant used to compute weights.
grid_ratio: tuple of int
Tuple of integers representing the ratio between the coarse and the fine grid.
dtype: np.dtype
Force to cast dtype for all matrices (interpolator.M may contain rationals).
Attributes
----------
dim: int
Number of dimensions to interpolate (same as interpolator.dim).
ghosts: tuple of int
Number of required ghosts.
n: tuple of int
Corresponds to 2*(ghosts+1), the number of required nodes to generate the
polynomial coefficients (same as interpolator.n).
N: int
Total number of input nodes N including ghosts (same as interpolator.N).
N = n0*n1*...*n[dim-1]
s: tuple of int
Corresponds to grid_ratio + 1, number of points of the subgrid in each directions.
Example for a grid ratio=(3,3), we have s=(4,4):
O=coarse grid nodes, X=fine grid nodes
Coarse grid: Fine grid:
^ O-----O ^ O X X O
| | | | X X X X
1 | | | 4 | X X X X
v O-----O v O X X O
<-----> <----->
1 4
S: int
Represents the number of fine grid points contained in a coarse grid cell.
S = s0*s1*...*s[dim-1]
gr: tuple of int
Corresponds to grid_ratio, number of points of the subgrid in each directions,
minus one.
Example for a grid ratio=(3,3), we have gr=(3,3) and s=(4,4):
O=coarse grid nodes, X=fine grid nodes, -=excluded find grid nodes
Coarse grid: Fine grid:
^ O-----O ^ O X X O ^
| | | gr0 | X X X - | s0
1 | | | v X X X - |
v O-----O O - - O v
<-----> <--->
1 gr1
GR: int
Represents the number of fine grid points contained in a coarse grid cell exluding
right most points.
GR = gr0*gr1*...*gr[dim-1]
W: np.ndarray
Pre computed weights to interpolate directly from coarse to fine grid.
Let F be the vector of N known coarse grid node values (including required
ghosts).
Let G be the vector of S unknown fine grid node values.
N
<--------->
X X X X X X ^ |g0| ^
X X X X X X | |f0| ^ |g1| |
X X X X X X | |f1| | |g2| |
X X X X X X | |f2| | |g3| |
W = X X X X X X | S F= |f3| | N G = W*F = |g4| | S
X X X X X X | |f4| | |g5| |
X X X X X X | |f5| v |g6| |
X X X X X X | |g7| |
X X X X X X v |g8| v
Will contain the same data type as intepolator.M if dtype is not passed,
else W will be computed from user given dtype.
Wr: np.ndarray
Reduced W that exludes rightmost output points of ndimensional output vector.
Pre computed weights to interpolate directly from coarse to inner fine grid.
Let F be the vector of N known coarse grid node values (including required ghosts).
Let G be the vector of GR unknown fine inner grid node values (see gr attribute).
N
<--------->
X X X X X X ^ |g0| ^
X X X X X X | |f0| ^ |g1| |
X X X X X X | |f1| | |g2| |
X X X X X X | |f2| | |g3| |
Wr = X X X X X X | GR F= |f3| | N G = W*F = |g4| | GR
X X X X X X | |f4| | |g5| |
X X X X X X | |f5| v |g6| |
X X X X X X | |g7| |
X X X X X X v |g8| v
Same data type as W.
"""
check_instance(interpolator, PolynomialInterpolator)
check_instance(grid_ratio, tuple, values=int, size=interpolator.dim, minval=1)
p = interpolator.p
n = interpolator.n
P = interpolator.P
N = interpolator.N
M = interpolator.M
ghosts = interpolator.ghosts
gr = grid_ratio
GR = np.prod(gr, dtype=np.int32)
del grid_ratio
s = tuple(gri + 1 for gri in gr)
S = np.prod(s, dtype=np.int32)
dim = interpolator.dim
key = ("PolynomialSubgridInterpolator", interpolator.key, gr, str(dtype))
self.p = p
self.P = p
self.s = s
self.S = S
self.n = n
self.N = N
self.gr = gr
self.GR = GR
self.ghosts = interpolator.ghosts
self.dim = dim
self.interpolator = interpolator
self.key = key
cache_file = interpolator.cache_file()
try:
data = load_data_from_cache(cache_file, key)
if data is not None:
W = data
if _check_matrix(W):
self.W = data
return
except Exception as e:
msg = f"Failed to load data from cache because:\n{e}"
warnings.warn(msg, HysopCacheWarning)
X = tuple(
np.asarray(tuple(sm.Rational(j, gr) for j in range(0, si)))
for i, (gr, si) in enumerate(zip(gr, s))
)
V = np.vander(X[0], N=p[0], increasing=True)
for i in range(1, dim):
Vi = np.vander(X[i], N=p[i], increasing=True)
V = np.multiply.outer(V, Vi)
even_axes = tuple(range(0, V.ndim, 2))
odd_axes = tuple(range(1, V.ndim, 2))
axes = even_axes + odd_axes
V = np.transpose(V, axes=axes).copy()
assert V.shape[:dim] == s
assert V.shape[dim:] == p
V = V.reshape((S, P))
W = V.dot(M)
assert W.shape == (S, N)
if dtype is not None:
W = W.astype(dtype)
update_cache(cache_file, key, W)
self.W = W
def __call__(self, F):
return self.W.dot(F.ravel()).reshape(self.s)
[docs]
def generate_subgrid_restrictor(self):
return PolynomialSubgridRestrictor(subgrid_interpolator=self)
@property
def Wr(self):
assert self.W.shape == (self.S, self.N)
view = (slice(0, -1, None),) * self.dim + (slice(None, None, None),) * self.dim
Wr = self.W.reshape(self.s + self.n)[view].reshape(self.GR, self.N).copy()
return Wr
def __hash__(self):
objs = (self.interpolator, self.gr)
return hash(objs)
[docs]
class PolynomialSubgridRestrictor:
def __init__(self, subgrid_interpolator):
"""
Create a PolynomialSubgridRestrictor from a PolynomialSubgridInterpolator.
Parameters
----------
subgrid_interpolator: PolynomialSubgridInterpolator
Interpolant used to compute restrictor weights.
Attributes
----------
g: tuple of int
Corresponds to (n+1)*s
G: int
Corresponds to g[0]*g[1]*...g[dim-1]
R: np.ndarray
Restrictor weights.
origin: tuple of int
Origin of the generated stencil.
Rr: np.ndarray
Restrictor weights excluding leftmost and rightmost points.
ghosts:
Corresponds to origin - 1, which is also Rr origin.
"""
check_instance(subgrid_interpolator, PolynomialSubgridInterpolator)
dim = subgrid_interpolator.dim
n = subgrid_interpolator.n
s = subgrid_interpolator.s
gr = subgrid_interpolator.gr
GR = subgrid_interpolator.GR
W = subgrid_interpolator.W
g = tuple(ni * gri + 1 for (ni, gri) in zip(n, gr))
G = np.prod(g, dtype=np.int64)
assert all(gi % 2 == 1 for gi in g)
origin = tuple(gi // 2 for gi in g)
gvals, gvars = tensor_symbol("g", g, origin)
I = 0
for idx in np.ndindex(*gr):
mask = tuple(
slice(i, max(2, i + ni * gri), gri) for (i, ni, gri) in zip(idx, n, gr)
)
target = tuple(gri - i for (i, gri) in zip(idx, gr))
F = gvals[mask]
Ii = W.dot(F.ravel()).reshape(s)[target]
I += Ii
R = np.ndarray(shape=g, dtype=object)
for idx in np.ndindex(*g):
R[idx] = I.coeff(gvals[idx])
view = (slice(1, -1, None),) * dim
Rr = R[view]
ghosts = tuple(oi - 1 for oi in origin)
self.g = g
self.G = G
self.R = R
self.Rr = Rr
self.origin = origin
self.ghosts = ghosts
self.n = n
self.s = s
self.GR = GR
self.gr = gr
self.subgrid_interpolator = subgrid_interpolator
if __name__ == "__main__":
np.set_printoptions(
precision=4,
linewidth=1e8,
threshold=1e8,
formatter={"float": lambda x: f"{x:+0.3f}"},
)
# 2D tests
grid_ratio = (2, 2)
F = [[1, 1], [1, 1]]
F = np.asarray(F)
print("Solving bilinear...")
PI = PolynomialInterpolator(dim=2, deg=1, fd=2, verbose=False)
GI0 = PI.generate_subgrid_interpolator(grid_ratio=grid_ratio)
GI1 = PI.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Bilinear (Rational)")
print(GI0(F))
print()
print("Bilinear (np.float64)")
print(GI1(F))
print()
F = [[0, 0, 0, 0], [0, 1, 1, 0], [0, 1, 1, 0], [0, 0, 0, 0]]
F = np.asarray(F)
print("Solving bicubic2...")
PI0 = PolynomialInterpolator(dim=2, deg=3, fd=2, verbose=False)
GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Bicubic (FDC2)")
print(GI0(F))
print()
print("Solving biquintic2...")
PI1 = PolynomialInterpolator(dim=2, deg=5, fd=2, verbose=False)
GI1 = PI1.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Biquintic (FDC2)")
print(GI1(F))
print()
F = [[0, 1, 1, 0], [0, 1, 1, 0]]
F = np.asarray(F)
print("Solving linear/cubic...")
PI0 = PolynomialInterpolator(dim=2, deg=(1, 3), fd=2, verbose=False)
GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Linear/Cubic (FDC2)")
print(GI0(F))
print()
F = [[0, 0], [1, 1], [1, 1], [0, 0]]
F = np.asarray(F)
print("Solving cubic/linear...")
PI0 = PolynomialInterpolator(dim=2, deg=(3, 1), fd=2, verbose=False)
GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Cubic/Linear (FDC2)")
print(GI0(F))
print()
F = [
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 1, 1, 0, 0],
[0, 0, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
]
F = np.asarray(F)
print("Solving bicubic4...")
PI0 = PolynomialInterpolator(dim=2, deg=3, fd=4, verbose=False)
GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Bicubic (FDC4)")
print(GI0(F))
print()
print("Solving biquintic4...")
PI1 = PolynomialInterpolator(dim=2, deg=5, fd=4, verbose=False)
GI1 = PI1.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Biquintic (FDC4)")
print(GI1(F))
print()
print("Solving biseptic2...")
PI2 = PolynomialInterpolator(
dim=2, deg=7, fd=2, verbose=False, approximative=(not has_flint)
)
GI2 = PI2.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Biseptic (FDC2)")
print(GI2(F))
print()
print("Solving binonic2...")
PI3 = PolynomialInterpolator(
dim=2, deg=9, fd=2, verbose=False, approximative=(not has_flint)
)
GI3 = PI3.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("Binonic (FDC2)")
print(GI3(F))
print()
print("Solving septic2/nonic2 ...")
PI4 = PolynomialInterpolator(
dim=2, deg=(7, 9), fd=2, verbose=False, approximative=(not has_flint)
)
GI4 = PI4.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("septic/nonic (FDC2)")
print(GI4(F))
print()
print("Solving septic2/quintic4 ...")
PI5 = PolynomialInterpolator(
dim=2, deg=(7, 5), fd=(2, 4), verbose=False, approximative=(not has_flint)
)
GI5 = PI5.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print("septic/nonic (FDC2/FDC4)")
print(GI5(F))
print()
# 3D test
grid_ratio = (2, 2, 2)
print("Solving trilinear...")
PI = PolynomialInterpolator(dim=3, deg=1, fd=2, verbose=False)
GI0 = PI.generate_subgrid_interpolator(grid_ratio=grid_ratio)
print()
print("Solving tricubic2...")
PI0 = PolynomialInterpolator(dim=3, deg=3, fd=2, verbose=False)
GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print()
print("Solving triquintic2...")
PI0 = PolynomialInterpolator(
dim=3, deg=5, fd=2, verbose=False, approximative=(not has_flint)
)
GI0 = PI0.generate_subgrid_interpolator(grid_ratio=grid_ratio, dtype=np.float64)
print()